% This code computes the First-best allocation of a RBC model with 

clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

isOctave && warning('off', 'Octave:array-as-logical');

load ../toDynare

tauk = eco.tauk;


alpha = eco.alpha(1);
beta = eco.beta;
delta = eco.delta(1);
kappa = 0;
chi = eco.chi;
phi = eco.phi;
G = solution.G;
coefB = 0;

rho_u = 0.1;
u0 = .01*(1-rho_u/(solution.R))/solution.G; 



file = 'code_dynare_RA_rules.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';

GsY = solution.G/solution.Y;
chi = 1;

BsY = 0;

dif = 1;


r = 1/beta - 1;
FK = r;
KsL = ((FK+delta)/alpha)^(1/(alpha - 1));
FL = (1-alpha)*KsL^alpha;
w = FL;
Ltot =  (chi*w)^phi;
K = KsL*Ltot;
Y = K^alpha*Ltot^(1-alpha);
I = delta*K;
Ctot = Y*(1-GsY) -I; 
G = GsY*Y;
A = (Ctot +G - w*Ltot)/r; % check a=K
x = Ctot - (1/chi)*Ltot^(1+1/phi)/(1+1/phi);


if eco.sigma==1
    util = @(c) log(c);
else
    util = @(c) (c^(1-eco.sigma)-1)/(1-eco.sigma);
end

%%%%%%%%%%%%%%%% Welfare  %%%%%%%%
W = util(x);



%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     variables and params     %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['var\n'] ; fprintf(fid, str);

str = ['u r w FK FL K Ltot Ctot Y I It\n'] ; fprintf(fid, str);
str = ['ut Kt Ct wt rt Yt Lt G x W;\n'] ; fprintf(fid, str);


str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        initialisation        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(chi,formatSpec),';\n']; fprintf(fid, str);

%%
equa = 1; % Numbering equations
tol = 0;  % threshold for transition, to avoid to consider very small transitions

str = ['\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       model definition       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['model;\n\n']; fprintf(fid, str);

    
str = ['x = Ctot - (1/chi)*Ltot^(1+1/phi)/(1+1/phi); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['(x^-gamma) =  beta*(1+r(+1))*(x(+1)^-gamma); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
%%%%%%%% Resource contraint
str = ['G = Gss*(1+u); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['I = K - (1-delta)*K(-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['G  + Ctot + I = Y;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%str = ['a = B + K; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Production side
str = ['FL = (1 - alpha)*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Wage
str = ['Ltot =  (chi*w)^phi; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%% Taxes
str = ['r = FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['w = FL;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Welfare
if eco.sigma==1
    str = ['W = log(x);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
else
    str = ['W =  (x',num2str(h),'^(1-gamma)-1)/(1-gamma) ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
end;



%%%%%%%%%%%%%%%% log-deviation variables %%%%%%%%%%%%%%%%
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(Ctot,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
%str = ['Bt = 100*(B/',num2str(B,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(w,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - (1/beta - 1));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(Ltot,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(Y,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['It = 100*(I/',num2str(delta*K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n'] ; fprintf(fid, str);
%%%%%%%%%%%%%%%% end of the model part


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



str = ['initval;\n\n'] ; fprintf(fid, str);
    
str = ['Ltot     = ', num2str(Ltot,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = 1/beta - 1;\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(FK,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(K,formatSpec),';\n'] ; fprintf(fid, str);
str = ['x     = ', num2str(x,formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(Ctot,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['I  = delta*K;\n'] ; fprintf(fid, str);

str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['It    = 0;\n'] ; fprintf(fid, str);

str = ['G     = ', num2str(G,formatSpec),';\n'] ; fprintf(fid, str);
str = ['W     = ', num2str(W,formatSpec),';\n'] ; fprintf(fid, str);

str = ['end;\n\n'] ; fprintf(fid, str);
%%%%%%%%%%% end of the steady state part

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       Aggregate shock        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['shocks;\n var eps; stderr ',num2str(u0,formatSpec),';\n end;\n'] ; fprintf(fid, str);


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%         Stoch simul          %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['stoch_simul (order=1,irf=500) ut Yt Ct It Lt Kt W;'] ; fprintf(fid, str);

isOctave && fflush(fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        Launch Dynare         %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    command = ['dynare ',file, ' noclearall'];
    eval(command);

name = ['fig',num2str(rho_u),'FB.mat'];
save (name,'ut_eps','Yt_eps','It_eps','Lt_eps','Ct_eps','Kt_eps','W','W_eps','Ctot','Ltot')




